<!DOCTYPE html>
<HTML lang = "en">
<HEAD>
  <meta charset="UTF-8"/>
  <meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
  <title>Uncertainty Programming, Generalized Uncertainty Quantification</title>
  

  <script type="text/x-mathjax-config">
    MathJax.Hub.Config({
      tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]},
      TeX: { equationNumbers: { autoNumber: "AMS" } }
    });
  </script>

  <script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
  </script>

  
<style>
pre.hljl {
    border: 1px solid #ccc;
    margin: 5px;
    padding: 5px;
    overflow-x: auto;
    color: rgb(68,68,68); background-color: rgb(251,251,251); }
pre.hljl > span.hljl-t { }
pre.hljl > span.hljl-w { }
pre.hljl > span.hljl-e { }
pre.hljl > span.hljl-eB { }
pre.hljl > span.hljl-o { }
pre.hljl > span.hljl-k { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kc { color: rgb(59,151,46); font-style: italic; }
pre.hljl > span.hljl-kd { color: rgb(214,102,97); font-style: italic; }
pre.hljl > span.hljl-kn { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kp { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kr { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kt { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-n { }
pre.hljl > span.hljl-na { }
pre.hljl > span.hljl-nb { }
pre.hljl > span.hljl-nbp { }
pre.hljl > span.hljl-nc { }
pre.hljl > span.hljl-ncB { }
pre.hljl > span.hljl-nd { color: rgb(214,102,97); }
pre.hljl > span.hljl-ne { }
pre.hljl > span.hljl-neB { }
pre.hljl > span.hljl-nf { color: rgb(66,102,213); }
pre.hljl > span.hljl-nfm { color: rgb(66,102,213); }
pre.hljl > span.hljl-np { }
pre.hljl > span.hljl-nl { }
pre.hljl > span.hljl-nn { }
pre.hljl > span.hljl-no { }
pre.hljl > span.hljl-nt { }
pre.hljl > span.hljl-nv { }
pre.hljl > span.hljl-nvc { }
pre.hljl > span.hljl-nvg { }
pre.hljl > span.hljl-nvi { }
pre.hljl > span.hljl-nvm { }
pre.hljl > span.hljl-l { }
pre.hljl > span.hljl-ld { color: rgb(148,91,176); font-style: italic; }
pre.hljl > span.hljl-s { color: rgb(201,61,57); }
pre.hljl > span.hljl-sa { color: rgb(201,61,57); }
pre.hljl > span.hljl-sb { color: rgb(201,61,57); }
pre.hljl > span.hljl-sc { color: rgb(201,61,57); }
pre.hljl > span.hljl-sd { color: rgb(201,61,57); }
pre.hljl > span.hljl-sdB { color: rgb(201,61,57); }
pre.hljl > span.hljl-sdC { color: rgb(201,61,57); }
pre.hljl > span.hljl-se { color: rgb(59,151,46); }
pre.hljl > span.hljl-sh { color: rgb(201,61,57); }
pre.hljl > span.hljl-si { }
pre.hljl > span.hljl-so { color: rgb(201,61,57); }
pre.hljl > span.hljl-sr { color: rgb(201,61,57); }
pre.hljl > span.hljl-ss { color: rgb(201,61,57); }
pre.hljl > span.hljl-ssB { color: rgb(201,61,57); }
pre.hljl > span.hljl-nB { color: rgb(59,151,46); }
pre.hljl > span.hljl-nbB { color: rgb(59,151,46); }
pre.hljl > span.hljl-nfB { color: rgb(59,151,46); }
pre.hljl > span.hljl-nh { color: rgb(59,151,46); }
pre.hljl > span.hljl-ni { color: rgb(59,151,46); }
pre.hljl > span.hljl-nil { color: rgb(59,151,46); }
pre.hljl > span.hljl-noB { color: rgb(59,151,46); }
pre.hljl > span.hljl-oB { color: rgb(102,102,102); font-weight: bold; }
pre.hljl > span.hljl-ow { color: rgb(102,102,102); font-weight: bold; }
pre.hljl > span.hljl-p { }
pre.hljl > span.hljl-c { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-ch { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cm { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cp { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cpB { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cs { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-csB { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-g { }
pre.hljl > span.hljl-gd { }
pre.hljl > span.hljl-ge { }
pre.hljl > span.hljl-geB { }
pre.hljl > span.hljl-gh { }
pre.hljl > span.hljl-gi { }
pre.hljl > span.hljl-go { }
pre.hljl > span.hljl-gp { }
pre.hljl > span.hljl-gs { }
pre.hljl > span.hljl-gsB { }
pre.hljl > span.hljl-gt { }
</style>



  <style type="text/css">
  @font-face {
  font-style: normal;
  font-weight: 300;
}
@font-face {
  font-style: normal;
  font-weight: 400;
}
@font-face {
  font-style: normal;
  font-weight: 600;
}
html {
  font-family: sans-serif; /* 1 */
  -ms-text-size-adjust: 100%; /* 2 */
  -webkit-text-size-adjust: 100%; /* 2 */
}
body {
  margin: 0;
}
article,
aside,
details,
figcaption,
figure,
footer,
header,
hgroup,
main,
menu,
nav,
section,
summary {
  display: block;
}
audio,
canvas,
progress,
video {
  display: inline-block; /* 1 */
  vertical-align: baseline; /* 2 */
}
audio:not([controls]) {
  display: none;
  height: 0;
}
[hidden],
template {
  display: none;
}
a:active,
a:hover {
  outline: 0;
}
abbr[title] {
  border-bottom: 1px dotted;
}
b,
strong {
  font-weight: bold;
}
dfn {
  font-style: italic;
}
h1 {
  font-size: 2em;
  margin: 0.67em 0;
}
mark {
  background: #ff0;
  color: #000;
}
small {
  font-size: 80%;
}
sub,
sup {
  font-size: 75%;
  line-height: 0;
  position: relative;
  vertical-align: baseline;
}
sup {
  top: -0.5em;
}
sub {
  bottom: -0.25em;
}
img {
  border: 0;
}
svg:not(:root) {
  overflow: hidden;
}
figure {
  margin: 1em 40px;
}
hr {
  -moz-box-sizing: content-box;
  box-sizing: content-box;
  height: 0;
}
pre {
  overflow: auto;
}
code,
kbd,
pre,
samp {
  font-family: monospace, monospace;
  font-size: 1em;
}
button,
input,
optgroup,
select,
textarea {
  color: inherit; /* 1 */
  font: inherit; /* 2 */
  margin: 0; /* 3 */
}
button {
  overflow: visible;
}
button,
select {
  text-transform: none;
}
button,
html input[type="button"], /* 1 */
input[type="reset"],
input[type="submit"] {
  -webkit-appearance: button; /* 2 */
  cursor: pointer; /* 3 */
}
button[disabled],
html input[disabled] {
  cursor: default;
}
button::-moz-focus-inner,
input::-moz-focus-inner {
  border: 0;
  padding: 0;
}
input {
  line-height: normal;
}
input[type="checkbox"],
input[type="radio"] {
  box-sizing: border-box; /* 1 */
  padding: 0; /* 2 */
}
input[type="number"]::-webkit-inner-spin-button,
input[type="number"]::-webkit-outer-spin-button {
  height: auto;
}
input[type="search"] {
  -webkit-appearance: textfield; /* 1 */
  -moz-box-sizing: content-box;
  -webkit-box-sizing: content-box; /* 2 */
  box-sizing: content-box;
}
input[type="search"]::-webkit-search-cancel-button,
input[type="search"]::-webkit-search-decoration {
  -webkit-appearance: none;
}
fieldset {
  border: 1px solid #c0c0c0;
  margin: 0 2px;
  padding: 0.35em 0.625em 0.75em;
}
legend {
  border: 0; /* 1 */
  padding: 0; /* 2 */
}
textarea {
  overflow: auto;
}
optgroup {
  font-weight: bold;
}
table {
  font-family: monospace, monospace;
  font-size : 0.8em;
  border-collapse: collapse;
  border-spacing: 0;
}
td,
th {
  padding: 0;
}
thead th {
    border-bottom: 1px solid black;
    background-color: white;
}
tr:nth-child(odd){
  background-color: rgb(248,248,248);
}


/*
* Skeleton V2.0.4
* Copyright 2014, Dave Gamache
* www.getskeleton.com
* Free to use under the MIT license.
* http://www.opensource.org/licenses/mit-license.php
* 12/29/2014
*/
.container {
  position: relative;
  width: 100%;
  max-width: 960px;
  margin: 0 auto;
  padding: 0 20px;
  box-sizing: border-box; }
.column,
.columns {
  width: 100%;
  float: left;
  box-sizing: border-box; }
@media (min-width: 400px) {
  .container {
    width: 85%;
    padding: 0; }
}
@media (min-width: 550px) {
  .container {
    width: 80%; }
  .column,
  .columns {
    margin-left: 4%; }
  .column:first-child,
  .columns:first-child {
    margin-left: 0; }

  .one.column,
  .one.columns                    { width: 4.66666666667%; }
  .two.columns                    { width: 13.3333333333%; }
  .three.columns                  { width: 22%;            }
  .four.columns                   { width: 30.6666666667%; }
  .five.columns                   { width: 39.3333333333%; }
  .six.columns                    { width: 48%;            }
  .seven.columns                  { width: 56.6666666667%; }
  .eight.columns                  { width: 65.3333333333%; }
  .nine.columns                   { width: 74.0%;          }
  .ten.columns                    { width: 82.6666666667%; }
  .eleven.columns                 { width: 91.3333333333%; }
  .twelve.columns                 { width: 100%; margin-left: 0; }

  .one-third.column               { width: 30.6666666667%; }
  .two-thirds.column              { width: 65.3333333333%; }

  .one-half.column                { width: 48%; }

  /* Offsets */
  .offset-by-one.column,
  .offset-by-one.columns          { margin-left: 8.66666666667%; }
  .offset-by-two.column,
  .offset-by-two.columns          { margin-left: 17.3333333333%; }
  .offset-by-three.column,
  .offset-by-three.columns        { margin-left: 26%;            }
  .offset-by-four.column,
  .offset-by-four.columns         { margin-left: 34.6666666667%; }
  .offset-by-five.column,
  .offset-by-five.columns         { margin-left: 43.3333333333%; }
  .offset-by-six.column,
  .offset-by-six.columns          { margin-left: 52%;            }
  .offset-by-seven.column,
  .offset-by-seven.columns        { margin-left: 60.6666666667%; }
  .offset-by-eight.column,
  .offset-by-eight.columns        { margin-left: 69.3333333333%; }
  .offset-by-nine.column,
  .offset-by-nine.columns         { margin-left: 78.0%;          }
  .offset-by-ten.column,
  .offset-by-ten.columns          { margin-left: 86.6666666667%; }
  .offset-by-eleven.column,
  .offset-by-eleven.columns       { margin-left: 95.3333333333%; }

  .offset-by-one-third.column,
  .offset-by-one-third.columns    { margin-left: 34.6666666667%; }
  .offset-by-two-thirds.column,
  .offset-by-two-thirds.columns   { margin-left: 69.3333333333%; }

  .offset-by-one-half.column,
  .offset-by-one-half.columns     { margin-left: 52%; }

}
html {
  font-size: 62.5%; }
body {
  font-size: 1.5em; /* currently ems cause chrome bug misinterpreting rems on body element */
  line-height: 1.6;
  font-weight: 400;
  font-family: "Raleway", "HelveticaNeue", "Helvetica Neue", Helvetica, Arial, sans-serif;
  color: #222; }
h1, h2, h3, h4, h5, h6 {
  margin-top: 0;
  margin-bottom: 2rem;
  font-weight: 300; }
h1 { font-size: 3.6rem; line-height: 1.2;  letter-spacing: -.1rem;}
h2 { font-size: 3.4rem; line-height: 1.25; letter-spacing: -.1rem; }
h3 { font-size: 3.2rem; line-height: 1.3;  letter-spacing: -.1rem; }
h4 { font-size: 2.8rem; line-height: 1.35; letter-spacing: -.08rem; }
h5 { font-size: 2.4rem; line-height: 1.5;  letter-spacing: -.05rem; }
h6 { font-size: 1.5rem; line-height: 1.6;  letter-spacing: 0; }

p {
  margin-top: 0; }
a {
  color: #1EAEDB; }
a:hover {
  color: #0FA0CE; }
.button,
button,
input[type="submit"],
input[type="reset"],
input[type="button"] {
  display: inline-block;
  height: 38px;
  padding: 0 30px;
  color: #555;
  text-align: center;
  font-size: 11px;
  font-weight: 600;
  line-height: 38px;
  letter-spacing: .1rem;
  text-transform: uppercase;
  text-decoration: none;
  white-space: nowrap;
  background-color: transparent;
  border-radius: 4px;
  border: 1px solid #bbb;
  cursor: pointer;
  box-sizing: border-box; }
.button:hover,
button:hover,
input[type="submit"]:hover,
input[type="reset"]:hover,
input[type="button"]:hover,
.button:focus,
button:focus,
input[type="submit"]:focus,
input[type="reset"]:focus,
input[type="button"]:focus {
  color: #333;
  border-color: #888;
  outline: 0; }
.button.button-primary,
button.button-primary,
input[type="submit"].button-primary,
input[type="reset"].button-primary,
input[type="button"].button-primary {
  color: #FFF;
  background-color: #33C3F0;
  border-color: #33C3F0; }
.button.button-primary:hover,
button.button-primary:hover,
input[type="submit"].button-primary:hover,
input[type="reset"].button-primary:hover,
input[type="button"].button-primary:hover,
.button.button-primary:focus,
button.button-primary:focus,
input[type="submit"].button-primary:focus,
input[type="reset"].button-primary:focus,
input[type="button"].button-primary:focus {
  color: #FFF;
  background-color: #1EAEDB;
  border-color: #1EAEDB; }
input[type="email"],
input[type="number"],
input[type="search"],
input[type="text"],
input[type="tel"],
input[type="url"],
input[type="password"],
textarea,
select {
  height: 38px;
  padding: 6px 10px; /* The 6px vertically centers text on FF, ignored by Webkit */
  background-color: #fff;
  border: 1px solid #D1D1D1;
  border-radius: 4px;
  box-shadow: none;
  box-sizing: border-box; }
/* Removes awkward default styles on some inputs for iOS */
input[type="email"],
input[type="number"],
input[type="search"],
input[type="text"],
input[type="tel"],
input[type="url"],
input[type="password"],
textarea {
  -webkit-appearance: none;
     -moz-appearance: none;
          appearance: none; }
textarea {
  min-height: 65px;
  padding-top: 6px;
  padding-bottom: 6px; }
input[type="email"]:focus,
input[type="number"]:focus,
input[type="search"]:focus,
input[type="text"]:focus,
input[type="tel"]:focus,
input[type="url"]:focus,
input[type="password"]:focus,
textarea:focus,
select:focus {
  border: 1px solid #33C3F0;
  outline: 0; }
label,
legend {
  display: block;
  margin-bottom: .5rem;
  font-weight: 600; }
fieldset {
  padding: 0;
  border-width: 0; }
input[type="checkbox"],
input[type="radio"] {
  display: inline; }
label > .label-body {
  display: inline-block;
  margin-left: .5rem;
  font-weight: normal; }
ul {
  list-style: circle; }
ol {
  list-style: decimal; }
ul ul,
ul ol,
ol ol,
ol ul {
  margin: 1.5rem 0 1.5rem 3rem;
  font-size: 90%; }
li > p {margin : 0;}
th,
td {
  padding: 12px 15px;
  text-align: left;
  border-bottom: 1px solid #E1E1E1; }
th:first-child,
td:first-child {
  padding-left: 0; }
th:last-child,
td:last-child {
  padding-right: 0; }
button,
.button {
  margin-bottom: 1rem; }
input,
textarea,
select,
fieldset {
  margin-bottom: 1.5rem; }
pre,
blockquote,
dl,
figure,
table,
p,
ul,
ol,
form {
  margin-bottom: 1.0rem; }
.u-full-width {
  width: 100%;
  box-sizing: border-box; }
.u-max-full-width {
  max-width: 100%;
  box-sizing: border-box; }
.u-pull-right {
  float: right; }
.u-pull-left {
  float: left; }
hr {
  margin-top: 3rem;
  margin-bottom: 3.5rem;
  border-width: 0;
  border-top: 1px solid #E1E1E1; }
.container:after,
.row:after,
.u-cf {
  content: "";
  display: table;
  clear: both; }

pre {
  display: block;
  padding: 9.5px;
  margin: 0 0 10px;
  font-size: 13px;
  line-height: 1.42857143;
  word-break: break-all;
  word-wrap: break-word;
  border: 1px solid #ccc;
  border-radius: 4px;
}

pre.hljl {
  margin: 0 0 10px;
  display: block;
  background: #f5f5f5;
  border-radius: 4px;
  padding : 5px;
}

pre.output {
  background: #ffffff;
}

pre.code {
  background: #ffffff;
}

pre.julia-error {
  color : red
}

code,
kbd,
pre,
samp {
  font-family: Menlo, Monaco, Consolas, "Courier New", monospace;
  font-size: 0.9em;
}


@media (min-width: 400px) {}
@media (min-width: 550px) {}
@media (min-width: 750px) {}
@media (min-width: 1000px) {}
@media (min-width: 1200px) {}

h1.title {margin-top : 20px}
img {max-width : 100%}
div.title {text-align: center;}

  </style>
</HEAD>

<BODY>
  <div class ="container">
    <div class = "row">
      <div class = "col-md-12 twelve columns">
        <div class="title">
          <h1 class="title">Uncertainty Programming, Generalized Uncertainty Quantification</h1>
          <h5>Chris Rackauckas</h5>
          <h5>December 12th, 2020</h5>
        </div>

        <p>In this lecture we will mix two separate topics: uncertainty quantification and adaptivity of algorithms. Using compiler-based tooling, similar to how automatic differentiation and probabilsitic programming toolchains, we will show how one can begin to pushforward uncertainties of a model or calculation. This leads to an idea of <em>uncertainty programming</em>, a term which is not in use but should be justified by these notes.</p>
<h2>What is Uncertainty Quantification?</h2>
<p>Uncertainty quantification is the identification and quantification of sources of uncertainty. In our training of a neural differential equation, we have seen that the question of uncertainty can quickly become muddled. Results are inexact because of:</p>
<ul>
<li><p>Truncation errors in the ODE solve</p>
</li>
<li><p>Truncation errors in the adjoint ODE solve</p>
</li>
<li><p>Truncation errors in the interpolation calculation</p>
</li>
<li><p>Numerical errors in every dot product along the way &#40;&#33;&#41;</p>
</li>
<li><p>Numerical errors in matrix multiplication and linear solving &#40;latter when implicit&#41;</p>
</li>
<li><p>Numerical errors in backpropogation</p>
</li>
<li><p>Measurement errors in our fitting data</p>
</li>
<li><p>Randomness in the optimizer &#40;when stochastic, like ADAM&#41;</p>
</li>
<li><p>What is the error in the model specification / model form?</p>
</li>
</ul>
<p>&quot;How correct is my model?&quot; is thus a very involved question, since you&#39;d have to konw that every source of uncertainty is contained. In some cases we have rigorous mathematical results proving bounds. In other cases, we need to find empirical ways to quantify what&#39;s going on using our known bounds.</p>
<h2>Some High Level UQ Techniques</h2>
<p>Two high level UQ techniques fall out of methodologies we have recently discussed. If we fit a model <span class="math">$f$</span> to data, be it a neural network, a neural ODE, or some physical ODE model, we can fit it probabilistically using the Bayesian estimation or probabilistic programming tools previously described. With this form of fitting, one can ask the question &quot;what are the likely results from the model given these parameter distributions?&quot;, which can then be answered through Monte Carlo sampling.</p>
<p>Another form of high level UQ is global sensitivity analysis, which gives a measurement for how much the output is going to vary over a wide range and thus relates uncertainties in the input to uncertainties in the output.</p>
<h2>Pushforward Methods for Uncertainties</h2>
<p>Instead of relying on expensive Monte Carlo methods for the pushforward of an uncertainty, we can derive a more programmatic approach to uncertainty quantification through the use of uncertain number arithmetic.</p>
<p>To start, let&#39;s first revive the old physics way to doing simple uncertainty quantification. If you have two numbers, <span class="math">$x = a \pm b$</span>, one might remember the rules like,</p>
<p class="math">\[
\alpha + a = (\alpha + a) \pm b
\]</p>
<p class="math">\[
\alpha a = \alpha a \pm |\alpha| b
\]</p>
<p>Let&#39;s investigate this a bit more and see if we can develop an arithmetic, like dual numbers, to then propogate through whole programs. This idea comes from the arthmetic on normally distributed random variables. If we interpret <span class="math">$x \sim N(a,b)$</span>, i.e. a normally distributed random variable with mean <span class="math">$a$</span> and standard deviation <span class="math">$b$</span>, then the distributions follow that:</p>
<p class="math">\[
\alpha + a \sim N(\alpha + a,b)
\]</p>
<p class="math">\[
\alpha a \sim N(\alpha a,|\alpha|b)
\]</p>
<p>From here we can begin to expand to multiple variables. If <span class="math">$f = Ax$</span> where <span class="math">$x ~ N(\mu,\Sigma)$</span> is a multidimensional random variable, then</p>
<p class="math">\[
E[f] = A\mu
\]</p>
<p>and</p>
<p class="math">\[
V[f] = A \Sigma A^T
\]</p>
<p>Now take a nonlinear <span class="math">$f(x)$</span>. By a Taylor expansion we have that</p>
<p class="math">\[
f(x) = f_0 + Jx + \ldots
\]</p>
<p>i.e. the linear approximation is <span class="math">$f_0 + Jx$</span> where <span class="math">$f_0 = f(\mu)$</span> and <span class="math">$J$</span> is the Jacobian matrix. If we do a pushforward on the linear approximation, we receive</p>
<p class="math">\[
f(x) ~ N(f(\mu),J\Sigma J^T)
\]</p>
<p>which gives the rules for the pushforward on any possible function through the linearization. But the linearization is the same as the forward differencing ones, meaning that we can augment existing tooling for forward-mode automatic differentiation to perform pushforwards of uncertain quantities. A library which does this is <a href="https://github.com/JuliaPhysics/Measurements.jl">Measurements.jl</a>. Note that this library additionally tracks the correlations between each variable so that the second order terms are accurate.</p>
<h3>Measurements.jl in Practice: Measurements on DifferentialEquations</h3>
<p>Since DifferentialEquations.jl takes in arbitrary number types, we can have it recompile to do the arithemetic of uncertainty propogation. For example, the following solves the pendulum of arbitrary amplitude with respect to uncertain parameters and initial condtions:</p>
<p class="math">\[
\ddot{\theta} + \frac{g}{L} \sin(\theta) = 0
\]</p>


<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>OrdinaryDiffEq</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>Measurements</span><span class='hljl-t'>
</span><span class='hljl-n'>g</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nfB'>9.79</span><span class='hljl-t'> </span><span class='hljl-oB'>±</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.02</span><span class='hljl-p'>;</span><span class='hljl-t'> </span><span class='hljl-cs'># Gravitational constants</span><span class='hljl-t'>
</span><span class='hljl-n'>L</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nfB'>1.00</span><span class='hljl-t'> </span><span class='hljl-oB'>±</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.01</span><span class='hljl-p'>;</span><span class='hljl-t'> </span><span class='hljl-cs'># Length of the pendulum</span><span class='hljl-t'>

</span><span class='hljl-cs'>#Initial Conditions</span><span class='hljl-t'>
</span><span class='hljl-n'>u₀</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-ni'>0</span><span class='hljl-t'> </span><span class='hljl-oB'>±</span><span class='hljl-t'> </span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>π</span><span class='hljl-t'> </span><span class='hljl-oB'>/</span><span class='hljl-t'> </span><span class='hljl-ni'>3</span><span class='hljl-t'> </span><span class='hljl-oB'>±</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.02</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-cs'># Initial speed and initial angle</span><span class='hljl-t'>
</span><span class='hljl-n'>tspan</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>6.3</span><span class='hljl-p'>)</span><span class='hljl-t'>

</span><span class='hljl-cs'>#Define the problem</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>simplependulum</span><span class='hljl-p'>(</span><span class='hljl-n'>du</span><span class='hljl-p'>,</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>t</span><span class='hljl-p'>)</span><span class='hljl-t'>
    </span><span class='hljl-n'>θ</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'>
    </span><span class='hljl-n'>dθ</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'>
    </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>dθ</span><span class='hljl-t'>
    </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-p'>(</span><span class='hljl-n'>g</span><span class='hljl-oB'>/</span><span class='hljl-n'>L</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-oB'>*</span><span class='hljl-t'> </span><span class='hljl-nf'>sin</span><span class='hljl-p'>(</span><span class='hljl-n'>θ</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>

</span><span class='hljl-cs'>#Pass to solvers</span><span class='hljl-t'>
</span><span class='hljl-n'>prob</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>ODEProblem</span><span class='hljl-p'>(</span><span class='hljl-n'>simplependulum</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>u₀</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>tspan</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>sol</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>solve</span><span class='hljl-p'>(</span><span class='hljl-n'>prob</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nf'>Tsit5</span><span class='hljl-p'>(),</span><span class='hljl-t'> </span><span class='hljl-n'>reltol</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nfB'>1e-6</span><span class='hljl-p'>)</span><span class='hljl-t'>

</span><span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>Plots</span><span class='hljl-t'>
</span><span class='hljl-nf'>plot</span><span class='hljl-p'>(</span><span class='hljl-n'>sol</span><span class='hljl-p'>,</span><span class='hljl-n'>plotdensity</span><span class='hljl-oB'>=</span><span class='hljl-ni'>100</span><span class='hljl-p'>,</span><span class='hljl-n'>vars</span><span class='hljl-oB'>=</span><span class='hljl-ni'>2</span><span class='hljl-p'>)</span>
</pre>


<img src=""  />

<p>From here it is clear that as the pendulum goes forward, the uncertainty grows since the exact period is unclear. Notice another nice feature is on display here: the <a href="http://docs.juliaplots.org/latest/recipes/">Plots.jl Recipe System</a>. The plotting just worked because the recipes are a type-recursive system. The three steps were:</p>
<ol>
<li><p>The DifferentialEquations solution recipe transformed the ODE solution into an array of measurement variables</p>
</li>
<li><p>The Measurements recipe transformed the measurement variables into an array of floats along with a series of error bars</p>
</li>
<li><p>This array of floats was recognized as a native format, and thus the plot was made.</p>
</li>
</ol>
<p>Note that this idea of discretizing distributions and pushing them through a full calculation can be done with more accuracy by using things like orthogonal polynomial expsnsions. This is the <em>polynomial chaos expansion</em> approach which we will not cover, but there is a <a href="https://github.com/timueh/PolyChaos.jl">PolyChaos.jl package</a> one can explore.</p>
<h2>Quantifying Numerical Uncertainty with Intervals</h2>
<p>While Measurements gives a sense of uncertainty quantification for unknown inputs, a different form of uncertainty quantification is that for floating point uncertainties. For example, when you calculate <span class="math">$sin(2.3)$</span> on your computer, this has an error in the approximation, and what if we wanted to push these errors forward to get an interval which bounds the possible values given the numerical uncertainty? This is done via <em>interval arithemetic</em>.</p>
<p>For this we can use <a href="https://juliaintervals.github.io/IntervalArithmetic.jl/latest/">IntervalArithmetic.jl</a>. The idea of interval arithmetic is to work rigorously on sets of real numbers, i.e.</p>
<p class="math">\[
[a,b] = \{x\in\mathbb{R} : a\leq x \leq b\}
\]</p>
<p>We can construct an interval given the <code>interval</code> method:</p>


<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>IntervalArithmetic</span><span class='hljl-t'>
</span><span class='hljl-n'>x</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>interval</span><span class='hljl-p'>(</span><span class='hljl-nfB'>0.1</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.2</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
&#91;0.1, 0.200001&#93;
</pre>


<p>or using the shorthand</p>


<pre class='hljl'>
<span class='hljl-nfB'>0.1</span><span class='hljl-oB'>..</span><span class='hljl-nfB'>0.2</span>
</pre>


<pre class="output">
&#91;0.0999999, 0.200001&#93;
</pre>


<p>Here the operator catches the constants at compile time and notices that 0.1 and 0.2 cannot be exactly represented in floating point numbers, and thus on the left side it rounds down by one floating point number and on the right it rounds up by one floating point number to make sure the set rigorously contains the correct value.</p>
<p>From here, non-monotone functions can propogate intervals like:</p>


<pre class='hljl'>
<span class='hljl-p'>[</span><span class='hljl-n'>a</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>b</span><span class='hljl-p'>]</span><span class='hljl-oB'>^</span><span class='hljl-ni'>2</span><span class='hljl-t'> </span><span class='hljl-oB'>:=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-n'>a</span><span class='hljl-oB'>^</span><span class='hljl-ni'>2</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>b</span><span class='hljl-oB'>^</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'>  </span><span class='hljl-k'>if</span><span class='hljl-t'> </span><span class='hljl-ni'>0</span><span class='hljl-t'> </span><span class='hljl-oB'>&lt;</span><span class='hljl-t'> </span><span class='hljl-n'>a</span><span class='hljl-t'> </span><span class='hljl-oB'>&lt;</span><span class='hljl-t'> </span><span class='hljl-n'>b</span><span class='hljl-t'>
          </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nf'>max</span><span class='hljl-p'>(</span><span class='hljl-n'>a</span><span class='hljl-oB'>^</span><span class='hljl-ni'>2</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>b</span><span class='hljl-oB'>^</span><span class='hljl-ni'>2</span><span class='hljl-p'>)]</span><span class='hljl-t'>  </span><span class='hljl-k'>if</span><span class='hljl-t'> </span><span class='hljl-n'>a</span><span class='hljl-t'> </span><span class='hljl-oB'>&lt;</span><span class='hljl-t'> </span><span class='hljl-ni'>0</span><span class='hljl-t'> </span><span class='hljl-oB'>&lt;</span><span class='hljl-t'> </span><span class='hljl-n'>b</span><span class='hljl-t'>
          </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-n'>b</span><span class='hljl-oB'>^</span><span class='hljl-ni'>2</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>a</span><span class='hljl-oB'>^</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-k'>if</span><span class='hljl-t'> </span><span class='hljl-n'>a</span><span class='hljl-t'> </span><span class='hljl-oB'>&lt;</span><span class='hljl-t'> </span><span class='hljl-n'>b</span><span class='hljl-t'> </span><span class='hljl-oB'>&lt;</span><span class='hljl-t'> </span><span class='hljl-ni'>0</span>
</pre>


<p>The rules can get fairly complicated and may need to be derived for each individual elementary function, but, just like automatic differentiation, recursion can be performed to get to a bottom of primatives which is known to then propogate forward the intervals.</p>
<p>Because this form of uncertainty quantification is rigorous, we can prove theorem on it. For example, let&#39;s say we want to show that</p>


<pre class='hljl'>
<span class='hljl-nf'>h</span><span class='hljl-p'>(</span><span class='hljl-n'>x</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>x</span><span class='hljl-oB'>^</span><span class='hljl-ni'>2</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-ni'>2</span>
</pre>


<pre class="output">
h &#40;generic function with 1 method&#41;
</pre>


<p>has no roots in <span class="math">$[3,4]$</span>. Since these are rigorous bounds, it holds that</p>


<pre class='hljl'>
<span class='hljl-nf'>h</span><span class='hljl-p'>(</span><span class='hljl-nfB'>3..4</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
&#91;7, 14&#93;
</pre>


<p>is a rigorous bound on the possible values, and thus it must not have a root in this interval.</p>
<h3>Problems with Interval Arithmetic</h3>
<p>Interval arithmetic is nice... but it can have some issues. It&#39;s rigorous but it&#39;s also conservative, meaning that the intervals can be much larger than one would expect given the actual uncertainties seen in practice. One phonomena which causes this can be seen by looking at that pendulum:</p>


<pre class='hljl'>
<span class='hljl-n'>g</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nfB'>9.77</span><span class='hljl-oB'>..</span><span class='hljl-nfB'>9.81</span><span class='hljl-p'>;</span><span class='hljl-t'> </span><span class='hljl-cs'># Gravitational constants</span><span class='hljl-t'>
</span><span class='hljl-n'>L</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.99</span><span class='hljl-oB'>..</span><span class='hljl-nfB'>1.01</span><span class='hljl-p'>;</span><span class='hljl-t'> </span><span class='hljl-cs'># Length of the pendulum</span><span class='hljl-t'>

</span><span class='hljl-cs'>#Initial Conditions</span><span class='hljl-t'>
</span><span class='hljl-n'>u₀</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-nfB'>0..0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-p'>((</span><span class='hljl-n'>π</span><span class='hljl-t'> </span><span class='hljl-oB'>/</span><span class='hljl-t'> </span><span class='hljl-ni'>3</span><span class='hljl-p'>)</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>0.02</span><span class='hljl-p'>)</span><span class='hljl-oB'>..</span><span class='hljl-p'>((</span><span class='hljl-n'>π</span><span class='hljl-t'> </span><span class='hljl-oB'>/</span><span class='hljl-t'> </span><span class='hljl-ni'>3</span><span class='hljl-p'>)</span><span class='hljl-oB'>+</span><span class='hljl-nfB'>0.02</span><span class='hljl-p'>)]</span><span class='hljl-t'> </span><span class='hljl-cs'># Initial speed and initial angle</span><span class='hljl-t'>
</span><span class='hljl-n'>tspan</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>6.3</span><span class='hljl-p'>)</span><span class='hljl-t'>

</span><span class='hljl-cs'>#Define the problem</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>simplependulum</span><span class='hljl-p'>(</span><span class='hljl-n'>du</span><span class='hljl-p'>,</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>t</span><span class='hljl-p'>)</span><span class='hljl-t'>
    </span><span class='hljl-n'>θ</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'>
    </span><span class='hljl-n'>dθ</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'>
    </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>dθ</span><span class='hljl-t'>
    </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-p'>(</span><span class='hljl-n'>g</span><span class='hljl-oB'>/</span><span class='hljl-n'>L</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-oB'>*</span><span class='hljl-t'> </span><span class='hljl-nf'>sin</span><span class='hljl-p'>(</span><span class='hljl-n'>θ</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>

</span><span class='hljl-cs'>#Pass to solvers</span><span class='hljl-t'>
</span><span class='hljl-n'>prob</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>ODEProblem</span><span class='hljl-p'>(</span><span class='hljl-n'>simplependulum</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>u₀</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>tspan</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>sol</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>solve</span><span class='hljl-p'>(</span><span class='hljl-n'>prob</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nf'>Tsit5</span><span class='hljl-p'>(),</span><span class='hljl-t'> </span><span class='hljl-n'>adaptive</span><span class='hljl-oB'>=</span><span class='hljl-kc'>false</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>dt</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>0.1</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>reltol</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nfB'>1e-6</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
retcode: Success
Interpolation: specialized 4th order &quot;free&quot; interpolation
t: 64-element Array&#123;Float64,1&#125;:
 0.0
 0.1
 0.2
 0.30000000000000004
 0.4
 0.5
 0.6
 0.7
 0.7999999999999999
 0.8999999999999999
 ⋮
 5.4999999999999964
 5.599999999999996
 5.699999999999996
 5.799999999999995
 5.899999999999995
 5.999999999999995
 6.099999999999994
 6.199999999999994
 6.3
u: 64-element Array&#123;Array&#123;IntervalArithmetic.Interval&#123;Float64&#125;,1&#125;,1&#125;:
 &#91;&#91;0, 0&#93;, &#91;1.02719, 1.0672&#93;&#93;
 &#91;&#91;0.0768938, 0.129145&#93;, &#91;0.649387, 1.34336&#93;&#93;
 &#91;&#91;-0.65616, 1.04792&#93;, &#91;-4.56529, 6.2509&#93;&#93;
 &#91;&#91;-17.5149, 18.0655&#93;, &#91;-12.0408, 13.5935&#93;&#93;
 &#91;&#91;-43.3001, 44.0059&#93;, &#91;-19.5521, 21.1047&#93;&#93;
 &#91;&#91;-74.7788, 75.6399&#93;, &#91;-27.0633, 28.616&#93;&#93;
 &#91;&#91;-111.952, 112.968&#93;, &#91;-34.5745, 36.1272&#93;&#93;
 &#91;&#91;-154.818, 155.989&#93;, &#91;-42.0857, 43.6384&#93;&#93;
 &#91;&#91;-203.377, 204.704&#93;, &#91;-49.597, 51.1497&#93;&#93;
 &#91;&#91;-257.63, 259.113&#93;, &#91;-57.1082, 58.6609&#93;&#93;
 ⋮
 &#91;&#91;-8908.08, 8916.71&#93;, &#91;-402.625, 404.178&#93;&#93;
 &#91;&#91;-9229.93, 9238.71&#93;, &#91;-410.136, 411.689&#93;&#93;
 &#91;&#91;-9557.48, 9566.42&#93;, &#91;-417.648, 419.2&#93;&#93;
 &#91;&#91;-9890.72, 9899.81&#93;, &#91;-425.159, 426.712&#93;&#93;
 &#91;&#91;-10229.7, 10238.9&#93;, &#91;-432.67, 434.223&#93;&#93;
 &#91;&#91;-10574.3, 10583.7&#93;, &#91;-440.181, 441.734&#93;&#93;
 &#91;&#91;-10924.6, 10934.2&#93;, &#91;-447.693, 449.245&#93;&#93;
 &#91;&#91;-11280.7, 11290.4&#93;, &#91;-455.204, 456.756&#93;&#93;
 &#91;&#91;-11642.4, 11652.2&#93;, &#91;-462.715, 464.268&#93;&#93;
</pre>


<p>While we start out with reasonably small intervals, it turns out that every operation is calculating &quot;what is the largest I could be? What is the smallest I could be?&quot;. Comparing these extremes at every operation means that, yes, by the end, given the uncertainty in the period, the solution lies in the interval <span class="math">$[-11642.4,11652.2]$</span>, but that&#39;s not a particularly helpful estimate&#33; This demonstrates the exponential explosion of interval estimates.</p>
<p>But note that part of why it got so large is because we started with &quot;such large&quot; intervals. If we only used this to measure the uncertainty of the floating point arithmetic, then the intervals are much better contained:</p>


<pre class='hljl'>
<span class='hljl-n'>g</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nfB'>9.8</span><span class='hljl-oB'>..</span><span class='hljl-nfB'>9.8</span><span class='hljl-p'>;</span><span class='hljl-t'> </span><span class='hljl-cs'># Gravitational constants</span><span class='hljl-t'>
</span><span class='hljl-n'>L</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nfB'>1.0</span><span class='hljl-oB'>..</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>;</span><span class='hljl-t'> </span><span class='hljl-cs'># Length of the pendulum</span><span class='hljl-t'>

</span><span class='hljl-cs'>#Initial Conditions</span><span class='hljl-t'>
</span><span class='hljl-n'>u₀</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-nfB'>0..0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-n'>π</span><span class='hljl-t'> </span><span class='hljl-oB'>/</span><span class='hljl-t'> </span><span class='hljl-ni'>3</span><span class='hljl-p'>)</span><span class='hljl-oB'>..</span><span class='hljl-p'>(</span><span class='hljl-n'>π</span><span class='hljl-t'> </span><span class='hljl-oB'>/</span><span class='hljl-t'> </span><span class='hljl-ni'>3</span><span class='hljl-p'>)]</span><span class='hljl-t'> </span><span class='hljl-cs'># Initial speed and initial angle</span><span class='hljl-t'>
</span><span class='hljl-n'>tspan</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>6.3</span><span class='hljl-p'>)</span><span class='hljl-t'>

</span><span class='hljl-cs'>#Define the problem</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>simplependulum</span><span class='hljl-p'>(</span><span class='hljl-n'>du</span><span class='hljl-p'>,</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>t</span><span class='hljl-p'>)</span><span class='hljl-t'>
    </span><span class='hljl-n'>θ</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'>
    </span><span class='hljl-n'>dθ</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'>
    </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>dθ</span><span class='hljl-t'>
    </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-p'>(</span><span class='hljl-n'>g</span><span class='hljl-oB'>/</span><span class='hljl-n'>L</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-oB'>*</span><span class='hljl-t'> </span><span class='hljl-nf'>sin</span><span class='hljl-p'>(</span><span class='hljl-n'>θ</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>

</span><span class='hljl-cs'>#Pass to solvers</span><span class='hljl-t'>
</span><span class='hljl-n'>prob</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>ODEProblem</span><span class='hljl-p'>(</span><span class='hljl-n'>simplependulum</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>u₀</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>tspan</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>sol</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>solve</span><span class='hljl-p'>(</span><span class='hljl-n'>prob</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nf'>Vern9</span><span class='hljl-p'>(),</span><span class='hljl-t'> </span><span class='hljl-n'>adaptive</span><span class='hljl-oB'>=</span><span class='hljl-kc'>false</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>dt</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>0.001</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>reltol</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nfB'>1e-6</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
retcode: Success
Interpolation: specialized 9th order lazy interpolation
t: 6301-element Array&#123;Float64,1&#125;:
 0.0
 0.001
 0.002
 0.003
 0.004
 0.005
 0.006
 0.007
 0.008
 0.009000000000000001
 ⋮
 6.292000000000436
 6.293000000000436
 6.294000000000437
 6.295000000000437
 6.296000000000437
 6.297000000000438
 6.298000000000438
 6.299000000000438
 6.3
u: 6301-element Array&#123;Array&#123;IntervalArithmetic.Interval&#123;Float64&#125;,1&#125;,1&#125;:
 &#91;&#91;0, 0&#93;, &#91;1.04719, 1.0472&#93;&#93;
 &#91;&#91;0.00104719, 0.0010472&#93;, &#91;1.04719, 1.0472&#93;&#93;
 &#91;&#91;0.00209438, 0.00209439&#93;, &#91;1.04717, 1.04718&#93;&#93;
 &#91;&#91;0.00314154, 0.00314155&#93;, &#91;1.04715, 1.04716&#93;&#93;
 &#91;&#91;0.00418868, 0.00418869&#93;, &#91;1.04711, 1.04712&#93;&#93;
 &#91;&#91;0.00523577, 0.00523578&#93;, &#91;1.04706, 1.04707&#93;&#93;
 &#91;&#91;0.00628281, 0.00628282&#93;, &#91;1.04701, 1.04702&#93;&#93;
 &#91;&#91;0.00732979, 0.0073298&#93;, &#91;1.04694, 1.04695&#93;&#93;
 &#91;&#91;0.0083767, 0.00837671&#93;, &#91;1.04686, 1.04687&#93;&#93;
 &#91;&#91;0.00942353, 0.00942354&#93;, &#91;1.04678, 1.04679&#93;&#93;
 ⋮
 &#91;&#91;-3.57381, 3.99756&#93;, &#91;-6.93095, 8.89265&#93;&#93;
 &#91;&#91;-3.58701, 4.01272&#93;, &#91;-6.94843, 8.91012&#93;&#93;
 &#91;&#91;-3.60024, 4.02791&#93;, &#91;-6.9659, 8.9276&#93;&#93;
 &#91;&#91;-3.6135, 4.04313&#93;, &#91;-6.98337, 8.94507&#93;&#93;
 &#91;&#91;-3.62679, 4.05838&#93;, &#91;-7.00085, 8.96254&#93;&#93;
 &#91;&#91;-3.64011, 4.07367&#93;, &#91;-7.01832, 8.98002&#93;&#93;
 &#91;&#91;-3.65346, 4.08898&#93;, &#91;-7.0358, 8.99749&#93;&#93;
 &#91;&#91;-3.66684, 4.10432&#93;, &#91;-7.05327, 9.01496&#93;&#93;
 &#91;&#91;-3.68026, 4.1197&#93;, &#91;-7.07074, 9.03244&#93;&#93;
</pre>


<h2>Contextual Uncertainty Quantification</h2>
<p>Those previous methods were non-contextual and worked directly through program modification. However, by not &quot;clumping&quot; interactions, uncertainty quantification can have overestimates like is seen with the interval growth. Thus, just like with reverse-mode AD, can we instead look for higher order uncertainty primatives on which to build such a system? When digging into reverse-mode AD, we saw that adjoint problems in engineering corresponded to the the reverse-mode rules for things like linear solve, eigenvalue problems, and the solution of ODEs. There does not seem to be a general analogue in the case of uncertainty quantification, but there is hope. Since this is a big enough field, people have found special cases where uncertainty can be quantified in interesting manners. Let&#39;s look specifically at ODEs.</p>
<h2>Quantifying Uncertainty in ODE Solves for Adaptivity</h2>
<p>However, in some sense, adaptive numerical methods work by embedding a form of uncertainty quantification. Let&#39;s take a look at the Bogaki-Shampine method for solving ODEs:</p>
<p>Notice that there&#39;s a <span class="math">$y_{n+1}$</span> and a <span class="math">$z_{n+1}$</span> for two separate solutions for the next time step. It so happens that <span class="math">$y$</span> is <span class="math">$\mathcal{O}(\Delta t^3)$</span> while <span class="math">$z$</span> is <span class="math">$\mathcal{O}(\Delta t^2)$</span>, meaning that <span class="math">$E = z_{n+1} - y_{n+1}$</span> is a <span class="math">$\mathcal{O}(\Delta t^2)$</span> estimate for the error in a given step, since the two must both be &quot;<span class="math">$\Delta t^2$</span> close enough&quot; to the true solution. Similarly, when we looked at the Dormand-Prince method in our homework, the tableau:</p>
<p><img src="https://user-images.githubusercontent.com/1814174/70629597-3b5cd380-1bf8-11ea-8a16-07bb5bdc0c3a.PNG" alt="DP tableau" /></p>
<p>had a second row as well, with the first being <span class="math">$\mathcal{O}(n^5)$</span> and the second being <span class="math">$\mathcal{O}(n^4)$</span>. Thus these Runge-Kutta methods naturally have error estimators. In standard usage, they are compared to the tolerances, like:</p>
<p class="math">\[
q = \frac{E}{\text{reltol}\max(z_n,z_{n+1}) + \text{abstol}}
\]</p>
<p>and when <span class="math">$q<1$</span>, the <span class="math">$\Delta t$</span> gives an error larger than the tolerances and so the step is rejected, decreased, and tried again. In many cases, one may control the error proportionally to this error estimator, i.e. the next <span class="math">$\Delta t$</span> is the product <span class="math">$q \Delta t$</span>.</p>
<p>That&#39;s all for adapting to a tolerance, but can we use this to propogate uncertainties? It turns out we can. This is known as the ProbInts method. Essentially, instead of an ODE, we can think of having solved a stochastic differential equation whose additive noise term of size which matches our error estimate. Specifically, adding a noise which is normally distributed with mean zero and standard deviation <span class="math">$(\Delta t)^{p}$</span>, where <span class="math">$p$</span> is the order of the adaptive error estimate &#40;i.e. the order of the lower approximation&#41;, is an approximation to the possible values that could have occured given the noise that was seen. By adding this at every step, we can then recover a distribution of possible solutions/trajectories.</p>
<h3>ProbInts in Action</h3>


<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>DiffEqUncertainty</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>fitz</span><span class='hljl-p'>(</span><span class='hljl-n'>du</span><span class='hljl-p'>,</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>t</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>V</span><span class='hljl-p'>,</span><span class='hljl-n'>R</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-t'>
  </span><span class='hljl-n'>a</span><span class='hljl-p'>,</span><span class='hljl-n'>b</span><span class='hljl-p'>,</span><span class='hljl-n'>c</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>p</span><span class='hljl-t'>
  </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>c</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>V</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>V</span><span class='hljl-oB'>^</span><span class='hljl-ni'>3</span><span class='hljl-oB'>/</span><span class='hljl-ni'>3</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-n'>R</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-p'>(</span><span class='hljl-ni'>1</span><span class='hljl-oB'>/</span><span class='hljl-n'>c</span><span class='hljl-p'>)</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>V</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'>  </span><span class='hljl-n'>a</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>b</span><span class='hljl-oB'>*</span><span class='hljl-n'>R</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-n'>u0</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>;</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>]</span><span class='hljl-t'>
</span><span class='hljl-n'>tspan</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>20.0</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>p</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-nfB'>0.2</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.2</span><span class='hljl-p'>,</span><span class='hljl-nfB'>3.0</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>prob</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>ODEProblem</span><span class='hljl-p'>(</span><span class='hljl-n'>fitz</span><span class='hljl-p'>,</span><span class='hljl-n'>u0</span><span class='hljl-p'>,</span><span class='hljl-n'>tspan</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>

</span><span class='hljl-n'>cb</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>AdaptiveProbIntsUncertainty</span><span class='hljl-p'>(</span><span class='hljl-ni'>5</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-cs'># 5th order method</span><span class='hljl-t'>
</span><span class='hljl-n'>sol</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>solve</span><span class='hljl-p'>(</span><span class='hljl-n'>prob</span><span class='hljl-p'>,</span><span class='hljl-nf'>Tsit5</span><span class='hljl-p'>())</span><span class='hljl-t'>
</span><span class='hljl-n'>ensemble_prob</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>EnsembleProblem</span><span class='hljl-p'>(</span><span class='hljl-n'>prob</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>sim</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>solve</span><span class='hljl-p'>(</span><span class='hljl-n'>ensemble_prob</span><span class='hljl-p'>,</span><span class='hljl-nf'>Tsit5</span><span class='hljl-p'>(),</span><span class='hljl-n'>trajectories</span><span class='hljl-oB'>=</span><span class='hljl-ni'>100</span><span class='hljl-p'>,</span><span class='hljl-n'>callback</span><span class='hljl-oB'>=</span><span class='hljl-n'>cb</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-nf'>plot</span><span class='hljl-p'>(</span><span class='hljl-n'>sim</span><span class='hljl-p'>,</span><span class='hljl-n'>vars</span><span class='hljl-oB'>=</span><span class='hljl-p'>(</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>1</span><span class='hljl-p'>),</span><span class='hljl-n'>linealpha</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>0.4</span><span class='hljl-p'>)</span>
</pre>


<img src=""  />


<pre class='hljl'>
<span class='hljl-n'>cb</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>AdaptiveProbIntsUncertainty</span><span class='hljl-p'>(</span><span class='hljl-ni'>5</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>sol</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>solve</span><span class='hljl-p'>(</span><span class='hljl-n'>prob</span><span class='hljl-p'>,</span><span class='hljl-nf'>Tsit5</span><span class='hljl-p'>())</span><span class='hljl-t'>
</span><span class='hljl-n'>ensemble_prob</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>EnsembleProblem</span><span class='hljl-p'>(</span><span class='hljl-n'>prob</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>sim</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>solve</span><span class='hljl-p'>(</span><span class='hljl-n'>ensemble_prob</span><span class='hljl-p'>,</span><span class='hljl-nf'>Tsit5</span><span class='hljl-p'>(),</span><span class='hljl-n'>trajectories</span><span class='hljl-oB'>=</span><span class='hljl-ni'>100</span><span class='hljl-p'>,</span><span class='hljl-n'>callback</span><span class='hljl-oB'>=</span><span class='hljl-n'>cb</span><span class='hljl-p'>,</span><span class='hljl-n'>abstol</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>1e-3</span><span class='hljl-p'>,</span><span class='hljl-n'>reltol</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>1e-1</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-nf'>plot</span><span class='hljl-p'>(</span><span class='hljl-n'>sim</span><span class='hljl-p'>,</span><span class='hljl-n'>vars</span><span class='hljl-oB'>=</span><span class='hljl-p'>(</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>1</span><span class='hljl-p'>),</span><span class='hljl-n'>linealpha</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>0.4</span><span class='hljl-p'>)</span>
</pre>


<img src=""  />

<p>Notice that while an interval estimate would have grown to allow all extremes together, this form keeps the trajectories alive, allowing them to fall back to the mode, which decreases the true uncertainty. This is thus a good explanation as to why general methods will overestimate uncertainty.</p>
<h2>Adjoints of Uncertainty and the Koopman Operator</h2>
<p>Everything that we&#39;ve demonstrated here so far can be thought of as &quot;forward mode uncertainty quantification&quot;. For every example we have constructed a method such that, for a known probability distribution in <code>x</code>, we build the probability distribution of the output of the program, and then compute quantites from that. On a dynamical system this pushforward of a measure is denoted by the Frobenius-Perron operator. With a pushforward operator <span class="math">$P$</span> and an initial uncertainty density <span class="math">$f$</span>, we can reprsent calculating the expected value of some cost function on the solution via:</p>
<p class="math">\[
\mathbb{E}[g(x)|X \sim Pf] = \int_{S(A)} P f(x) g(x) dx
\]</p>
<p>where <span class="math">$S$</span> is the program, i.e. <span class="math">$S(A)$</span> is the total set of points by pushing every value of <span class="math">$A$</span> through our program, and  <span class="math">$P f(x)$</span> is the pushforward operator applied to the probability distribution. What this means is that, to calculate the expectation on the output of our program, like to calculate the mean value of the ODE&#39;s solution given uncertainty in the parameters, we can pushforward the probability distribution to construct <span class="math">$Pf$</span> and on this probability distribution calculate the expected value of some <span class="math">$g$</span> cost function on the solution.</p>
<p>The problem, as seen earlier, is that pushing forward entire probability distributions is a fairly expensive process. We can instead think about doing the adjoint to this cost function, i.e. pulling back the cost function and computing it on the initial density. In terms of inner product notation, this would be doing:</p>
<p class="math">\[
\langle Pf,g \rangle = \langle f, Ug \rangle
\]</p>
<p>meaning <span class="math">$U$</span> is the adjoint operator to the pushforward <span class="math">$P$</span>. This operator is known as the Koopman operator. There are many properties one can use about the Koopman operator, one special property being it&#39;s a linear operator on the space of observables, but it also gives a nice expression for computing uncertainty expectations. Using the Koopman operator, we can rewrite the expectation as:</p>
<p class="math">\[
\mathbb{E}[g(x)|X \sim Pf] = \mathbb{E}[Ug(x)|X \sim f]
\]</p>
<p>or perform the integral on the pullback of the cost function, i.e.</p>
<p class="math">\[
\mathbb{E}[g(x)|X \sim f] = \int_A Ug(x) f(x) dx
\]</p>
<p>In images it looks like:</p>
<p><img src="https://user-images.githubusercontent.com/1814174/102001466-a7b55b80-3cc0-11eb-9208-0f751fdca590.PNG" alt="Koopman vs FP" /></p>
<p>This expression gives us a fast way to compute expectations on the program output without having to compute the full uncertainty distribution on the output. This can thus be used for <em>optimization under uncertainty</em>, i.e. the optimization of loss functions with respect to expectations of the program&#39;s output under the assumption of given input uncertainty distributions. For more information, see <a href="https://arxiv.org/abs/2008.08737">The Koopman Expectation: An Operator Theoretic Method for Efficient Analysis and Optimization of Uncertain Hybrid Dynamical Systems</a>.</p>


        <HR/>
        <div class="footer">
          <p>
            Published from <a href="uncertainty_programming.jmd">uncertainty_programming.jmd</a>
            using <a href="http://github.com/JunoLab/Weave.jl">Weave.jl</a> v0.10.6 on 2020-12-12.
          </p>
        </div>
      </div>
    </div>
  </div>
</BODY>

</HTML>
